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ABSTRACT 

We have accurately measured the bispectrum for four scale-free models of 
structure formation with the spectral index n — 1, 0, —1, and —2. The mea- 
surement is based on a new method that can effectively eliminate the alias and 
numerical artifacts, and reliably extend the analysis into the strongly non-linear 
regime. The work makes use of a set of state-of-the art N-body simulations 
that have significantly increased the resolution range compared with the previ- 
ous studies on the subject. With these measured results, we demonstrate that 
the measured bispectrum depends on the shape and size of fc-triangle even in the 
strongly nonlinear regime. It increases with wavenumber and decreases with the 
spectral index. These results are in contrast with the hypothesis that the reduced 
bispectrum is a constant in the strongly non-linear regime. We also show that 
the fitting formula of Scoccimarro & Frieman (1999) does not describe our simu- 
lation results well (with a typical error about 40 percent). In the end, we present 
a new fitting formula for the reduced bispectrum that is valid for —2 < n < 
with a typical error of 10 percent only. 

Subject headings: cosmology: theory - galaxies:clusters:general - large-scale struc- 
ture of universe -met hods: N-body simulations 

1. Introduction 

Large-scale structures in the Universe are thought to arise from small primordial fluc- 
tuations through gravitational amplification. It is known that gravitational clustering is a 
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non-linear process. When the density fluctuations are sufficiently small, the evolution of the 
structures can be studied using perturbation theory (PT). With the growth of the fluctu- 
ation, even for an initially Gaussian fluctuation, nonlinear gravitational instability induces 
non-Gaussian signatures in the density field. In the weakly non-linear regime, leading order 
(tree-level) perturbation theory (Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1994a; 
Bernardeau et al. 1994b; Lokas et al. 1995; Gaztahaga & Baugh 1995; Baugh, Gaztahaga & 
Efstathiou 1995; Bouchet et al. 1995) can describe the clustering properties successfully. As 
one approaches smaller scales, the loop corrections to the tree-level results are expected to be- 
come important (Scoccimarro & Frieman 1996; Scoccimarro 1997). In the non-linear regime, 
numerical simulations must be applied to follow the development of the cosmic structures. 

The n-point correlation functions have been widely used as a powerful tool for quantify- 
ing the statistical properties of a density field both in theoretical models and observational 
catalogs (Peebles 1980). For a Gaussian random field, the two-point correlation function 
(2PCF) or its Fourier transform, the power spectrum P(k) can completely characterize its 
statistical properties, with all higher-order (connected) correlation functions being zero. It 
requires the higher order correlation functions to describe the statistical properties of the 
non-Gaussian distribution resulting from gravitational instability (Peebles 1980; Fry 1984; 
Bernardeau et al. 2002 for an excellent review and references therein). 

The bispectrum, the three-point correlation function (3PCF) in Fourier space, is the 
lowest order statistic that probes the shape of large-scale structures generated by the gravi- 
tational clustering (Peebles 1980). Theoretical models of weakly non-linear 3PCF have been 
studied well in the literature based on PT. PT can describe properties of dark matter on large 
scales > 10 h^Mpc. It predicts that the 3PCF depends on the shape of the linear power 
spectrum and on the shape of the triangle configuration both in real space (Jing, Borner 
& Valdarnini 1995; Jing & Borner 1997; Frieman & Gaztahaga 1999; Barriga & Gaztahaga 
2002) and in Fourier space (Fry 1984; Scoccimarro et al. 1998; Scoccimarro et al. 1999). 
When the galaxy bias is considered, the bispectrum of galaxies contains information on the 
primordial fluctuation and on galaxy biasing (Fry 1994; Fry & Gaztahaga 1993; Hivon et 
al. 1995; Mo, Jing, & White 1997; Matarrese et al. 1997; Verde et al. 2002). Measuring 
the galaxy bispectrum on large scales can help break the degeneracy between the linear bias 
and the matter parameter Q m usually present in the dynamical analysis of galaxy redshift 
surveys (Fry 1984; Hivon et al. 1995; Matarrese et al. 1997; Verde et al. 1998; Scoccimarro 
et al. 1998). Several authors have started to measure the bias parameters from current large 
galaxy surveys (Frieman & Gaztahaga 1999 and Gaztahaga & Freiman 1994 for the APM 
galaxies; Scoccimarro et al. 2001c for IRAS galaxies; Verde et al. 2002 and Jing & Borner 
2003 for the 2dFGRS galaxies; Kayo et al. 2004 for SDSS galaxies). 
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A quantitative modeling of the 3PCF or bispectrum in the nonlinear regime is more 
challenging. There have been attempts to predict the 3PCF based on the so-called halo 
model (Ma & Fry 2000, Scoccimarro et al. 2001b, Takada & Jain 2003a, Wang et al. 2004). 
In their detailed modeling, Takada & Jain (2003a) found that the halo model prediction 
for the 3PCF agrees with the simulation result of Jing & Borner (1998) both in linear and 
strongly non-linear regimes, but fails on intermediate non-linear scales (~ 1 h~ 1 Mpc). They 
also pointed out that the 3PCF at the intermediate scales is very sensitive to the outer radial 
cut of halos, which could be the reason for the failure of the halo model. On the other hand, 
N-body simulations have been widely used to study the 3PCF or bispectrum in the nonlinear 
regime (Davis et al. 1985; Efstathiou et al. 1988). Based on extensive studies with N-body 
simulations, a fitting formula for the bispectrum was proposed for the scale-free models by 
Scoccimarro & Frieman(1999, hereafter SF99), and then extended for the cold dark matter 
(CDM) models by Scoccimarro & Couchman (2001). 

The fitting formula of Scoccimarro & Couchman (2001) was applied to calculating the 
skewness of the convergence field in the weak gravitational lensing survey (Van Waerbeke 
et al. 2001; Hamana et al. 2002). Nowadays weak lensing surveys have become detailed to 
start measuring the skewness of the lensing shear field (Pen et al. 2003; Zhang et al. 2003), 
and will soon become big enough to measure the three-point correlation function of cosmic 
shears (Schneider et al. 2003; Takada & Jain 2003b). It is important to reliably predict the 
skewness and the three-point correlation function of cosmic shears in cosmological models, 
as the observational determinations of these quantities are expected to yield a measure of 
the cosmological density parameter (Bartelmann & Schneider 2001 for an excellent review). 
Therefore, the motivation is high to derive an accurate model for the bispectrum from linear 
to non-linear scales, and present them in a form useful for the weak lensing survey analysis. 

In this paper, we present such a study to investigate how the non-linear evolution of the 
bispectrum proceeds with a set of scale-free simulations of 512 3 particles. The simulations 
have the initial spectral index n — 1, 0, —1, or —2. We find that the bispectrum depends 
not only on n but also on the shape and size of the triangle even in the strong non-linear 
regime. Our results show that the formula of SF99 cannot accurately describe the properties 
of the bispectrum in the non-linear regime. Comparing our measured non-linear bispectrum 
with the weakly non-linear bispectrum obtained from the second order PT (hereafter PT2), 
we have arrived at a new formula for the bispectrum that is significantly more accurate than 
that of SF99. 

In this analysis, we have carefully examined possible effects on the bispectrum mea- 
surement of the numerical artifacts, such as finite box size, force softening, and particle 
discreteness. We have also closely paid attention to the effect caused by the mass assign- 
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ment using the Fast Fourier Transform (FFT, i.e. the alias effect). In order to get the 
true power spectrum and bispectrum from the simulation, we have developed a procedure 
to correct for the numerical and alias effects. 

This paper is organized as follows. In section 2, we give a brief overview of the self- 
similar evolution of the density field. In section 3, we describe the numerical simulations 
used. In section 4, we introduce the power spectrum and the bispectrum, and outline the 
possible numerical artifacts in the measurements of the power spectrum and the bispectrum. 
In section 5, we describe our method to measure the bispectrum. We present the measured 
bispectrum, and give our fitting formula for the bispectrum in section 6. The last section 
contains our summary. 



2. Self-similarity 

Gravitational clustering from initial conditions has provided one of the classic problems 
in cosmology. To achieve a self-similar evolution, according to Peebles (1980) and Efstathiou 
et al.(1988), two conditions must be fulfilled: (I) The background cosmology should not 
contain any characteristic scales, thus the universe must be an Einstein-de Sitter one, and 
(2) The initial density field should have no characteristic length scale, thus its initial power 
spectrum must be a power law. 

For a self-similar clustering pattern, its physical properties remain the same during its 
evolution when the length scale R is scaled by a characteristic scale R . A simple choice of R 
for a self-similar clustering is the scale at which the fluctuations begin to become non-linear, 
i.e., the variance of the linear density field smoothed on this scale is unity, 

a 2 (R , a) = 1 . (I) 

The variance is defined by 

/dk 
Al(k,a)\W(kR)\ 2 -. (2) 

where W(x) is the Fourier transform of a window function (usually a top hat or Gaussian), 
A 2 (k) = do 2 jdXwk = (V fl /(27r) 3 )47rk 3 P(k) is the contribution to the fractional density vari- 
ance per unit In A; (V^ is a normalization volume), and a(t) is the expansion scale factor. For 
the scale-free initial power spectrum, P L oc a 2 k n , the characteristic scale satisfies 

R <xa(t)^. (3) 

The characteristic wavenumber k can be chosen to be Rq \ so k oc a(ty 2 /^ +n \ With the 
characteristic scale _R , all statistical measures of the density field can be expressed as a 
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similarity solution that is independent of time 

f(R,t) = g(R/Ro) or f(k,t) = g(kR Q ), (4) 

(Peebles 1980; Efstathiou et al. 1988; Colombi, Bouchet, & Hernquist 1996; Jain & Bertschinger 
1998). 

3. The numerical simulation 

We study the scale-free models that assume an Einstein-de Sitter universe (i.e. f2 = 1 
and A = 0), and a power- law P(k) oc k n with n being 1, 0, —1, and —2 respectively for 
the linear density power spectrum. For each model, we have one simulation of 512 3 particles 
that was produced by one of the authors (YPJ). The current simulations are constructed in 
a similar way as the scale-free simulation sample in Jing (1998), but have higher force and 
mass resolutions. They were generated with a parallel-vectorized P 3 M (i.e. Particle Particle 
Particle Mesh) code (Jing & Suto 2002) at the National Astronomical Observatory of Japan. 
The gravitational force is softened with the S2 form of Hockney and Eastwood (1981) with 
the softening parameter rj = 1 x 10~ 4 L (L is the simulation box size), so the force becomes 
Newtonian when the separation is larger than rj. The simulations are evolved for 2000 time 
steps with a total of ten (n=l, 0, —2) or eleven (n= —1) outputs at a constant logarithmic 
interval (A log a) in the scale factor a. Table 1 summarizes the parameters that are relevant 
to the discussion in the current work. 



4. Power and bispectrum 
4.1. Basic theory 

Let p(r) be the cosmic density field, with the mean density p. The density field can be 
represented by a dimensionless field 5(f) (which is usually referred to as the density contrast) 

m = ^ (5) 

P 

Based on the cosmological principle, we expect p(f) to be periodic in some large rectangular 
volume V^. Its Fourier transformation is then defined by 

S(k) = / 5(r)e ir ^dr. (6) 
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The density field of the simulation is periodic at the box size L. The requirement of 
periodicity restricts the allowed wavenumbers to harmonic boundary conditions 

k x = nk b , (k b = —,n = 2, -1, 0, 1, 2, 3 • • •), (7) 

with similar expressions for k y and k z . 

The power spectrum P(k) and the bispectrum B i23 = B(k±, k 2 , h) are defined as 

^(S(h)6(k 2 )) = S Dirac (kj + k 2 )P(k), 
(5(h) 5(k 2 ) 5(h)) = 5 Dirac (h + h + h)B 123 , (8) 

where (• • •) means ensemble average, Spirac is the Dirac delta, and the 5m rac (h + h + h) 
implies that the bispectrum is defined for configurations of wavenumbers that form closed 
triangles in /c-space. There are many ways to express the shape of a triangle. For a triangle 
with h, k 2 , h, and \h\ > \k 2 \, we can parameterize its shape by k, v, 9 as: 

k = \h\, v — -^-, 9 = arccos ( — ^r~\ . (9) 
\h\ \\h\\h\J 

In the tree-level PT, the bispectrum can be expressed as follows: 

B 123 = 2 F 2 (h, k 2 )P x P 2 + eye, (10) 
where Pi = P(h) (i = 1, 2, 3) , and F 2 (h, k 2 ) is the kernel function, 

T?(t t\ 5 , 1 ki-h f h k 2 \ 2 [h-k 2 \ 
For convenience we can define the reduced bispectrum Q as 

Q^M,k^ PiP2 + B p 2 + PsPl - < 12 > 

According to Eq. (9), Q can be expressed as a function of k, v, and 9. 



4.2. Measuring the bispectrum 

The Fourier modes of a particle distribution can be determined exactly using the ex- 
pression (Peebles 1980) 

1 N - 

i=i 
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Owing to the periodic boundary condition in the simulation, wavenumbers are restricted 
to the form defined by Eq. (7). It is inefficient to use Eq. (13) to compute the Fourier 
transformation for a simulation where both iV and the mode number considered are large. 
Here we use the Fast Fourier Transform (FFT) technique to compute S(k). In this case, 
there is an upper limit for k that is imposed by the finite sampling of the density field at the 
FFT mesh points, which is called the Nyquist wavenumber, 

k Ny = (14) 

where Ax = L/N m is the mesh spacing and N m is the dimension of the mesh. Then the 
power spectrum and bispectrum can be estimated through averaging all modes of which the 
wavenumbers ki in thin shells (fcj ± Ak, Ak <C fcj, i=l, 2, 3) satisfy Eq. (8): 

p (h) = - (tf(fciW*2))<W(fci + fc2), 

B 123 = ^ V (S(k 1 )5(k 2 )5(h))S Dl rac(ki + k 2 + h), (15) 

fcl,fe 2 ,fc36* 

where <£> is the set composed by all wavenumbers k\ (with k 2 = —hi) in the thin shell (ki±Ak, 
Ak <C hi), and ^ is the set composed by all triangles with the same shapes formed by k\ 
and k 2 (with ks = —k 2 — k\) in their thin shells respectively, m and m' are the numbers of 
the pairs and the triangles to be averaged. 

Although it is very efficient to compute S(k) with FFT, it is important to remember the 
numerical limitations caused by FFT. As shown in Jing (1992; 2004), the mass assignment 
onto a grid for FFT has two effects: the smoothing effect and the sampling effect. The 
smoothing effect has been considered by many authors in previous bistpectrum measurements 
(e.g. Scoccimarro et al. 1998, their Appendix), but the sampling effect has not. Because 
both effects are coupled, the smoothing effect cannot be fully corrected if the sampling effect 
is not considered. These effects must be taken into account for a precision analysis such as 
the current work. First, the particle distribution must be sampled at the FFT grid points. 
Here we adopt the Nearest-Grid-Point (NGP) (Efstathiou et al. 1985) approximation to 
assign the particle mass to the grid. This generally leads to a smoothing of the density 
field on the scale of ~ Ax in the coordinate space. Higher order mass assignment schemes, 
e.g. cloud-in-cell (CIC) or triangular shaped cloud (TSC), cannot avoid the smoothing issue 
either; in fact they increase the smoothing even to a larger scale. In the next subsection, we 
will show that the mass assignment with NGP has the smoothing effect for k > k Ny /3 which 
essentially limits the usable Fourier modes to k < k Ny /3 . For a 3D FFT with N m = 1024, k 
is thus limited to < 160/cfe that is not enough for fully exploring the non-linear properties. A 
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much bigger N m would require a huge computer resource for the FFT computation. In the 
next section, we will show that we can overcome this problem effectively with a 2D FFT. 



When measuring the power spectrum and bispectrum from a simulation we must take 
account for the numerical artifacts, such as the discreteness effect, the finite box size, and 
the force softening. These limit the dynamical range of the simulation, and thus affect the 
measured power spectrum and bispectrum. Other important effect is that introduced by the 
FFT. Here we will address how to correct and/or account for these numerical artifacts. 



Since the power spectrum and bispectrum are measured from simulations with a finite 
number of particles, we need to correct for the discreteness effect arising from the Poisson 
shot noise. According to Peebles (1980), we divide the volume into infinitesimal elements 
{<iVi} with rii objects inside dV{. The over-density can be written as: 5i = (rii — n)/n (n is 
the mean number density), and its Fourier transformation can be expressed as: 



where N = nV M , is the number of particles in V^. Since dVi is taken so small that rii is either 
or 1, we have rii = nf = n\ = ■ ■ ■ . 



4.3. Numerical effects on the power spectrum and bispectrum 



4-3.1. Discreteness effects 




(16) 




= (S(h)S*(k 2 )) + -S Dirac (k 1: k 2 ) . (17) 



Finally, we get the desired result: 



P(A;) = (|^)| 2 ) 



1 



(18) 



N 



iFrom Eq. (18) we know that the discreteness (or shot noise) effect gives an additional 
term 1/N to the power spectrum. We can correct for the shot noise in the power spectrum 
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easily. In analogy with the power spectrum, we write the bispectrum : 
(5 d (h)5 d (k 2 )5 d (k 3 )} = (5(h)5(k 2 )5(h)) 

+^[P(ki) + P(h) + P(h)} + ^- 2 , (19) 
The bispectrum with the shot noise removed can be expressed as: 
B(h,k 2 ,k 3 ) = (5 d (h)5 d (k 2 )5 d (k 3 )) 

-^[P(h) + P(k 2 ) + P(h)]-^- 2 - (20) 



4-3.2. Force softening and box size 

In the simulation, in order to suppress two-body encounters, a softening must be applied 
when calculating the gravitational interaction. This induces an error in the integration of 
particle trajectories at small scale. We must impose some constraint on the scale below 
which the numerical effect dominates the clustering in the simulation. The cutoff can be a 
few times of the softening length. 

The box size of the simulation is finite, so there is no clustering power beyond the 
simulation box. On one hand, this results in a limited number of Fourier modes at k > k^ 
which can influence the accuracy of measuring the clustering spectra. On the other hand, 
this large-scale cutoff may also affect the clustering on scale much smaller than the box size, 
because the coupling between different scales can be important on non-linear and quasilinear 
scales. 

Both the force softening and the box-size cutoff are expected to break down the scaling 
property of the self-similar evolution. We will use the expected scaling to quantify these 
numerical artifacts. 



4-3.3. Mass assignment 

When doing the FFT, we first need to collect density values on grids (usually called 
mass assignment). The mass assignment in fact is equivalent to convolving the density field 
by one chosen function W(r) and sampling the convolved density on a finite number of grid 
points, therefore the FFT of p(r g ) generally is not equal to the FT of p(r). In this work, 
we have adopted the NGP (the nearest grid point) scheme to assign particles to the mesh. 
The finite sampling of the convolved density results in the summation of the aliased power 
spectrum or bispectrum. 
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We can correct the alias effect on the power spectrum through a theoretical calculation 
(Jing 1992, 2004; Baugh & Efstathiou 1994; Smith et al. 2003). But it is complicated to cor- 
rect for the alias effect on the bispectrum through a theoretical calculation. Fortunately, we 
find that the 2D statistical properties are identical to the 3D statistical properties in Fourier 
space (§5). Thus we can make use of the 2D density field instead of the 3D density field 
when calculating the bispectrum at small scales, because the 2D FFT needs less computer 
memory than the 3D case. In fact, the larger number of grid points N m for the FFT, the 
smaller scale where the alias effect takes place. The 2D FFT can overcome the limitation of 
the computation and involve little alias effect. Thus we can extend the measurement of the 
bispectrum to very small scales with little numerical artifact. 



5. The method 

Let S^nif) be the over-density in 3D real space, and 52d(x) the over-density in 2D real 
space, 

Po 

5 2D (x) = -. (21) 

where p(r) and u(x) are the density field of the 3D real space and the 2D real space respec- 
tively, po and Uq are the mean density correspondingly. The 2D density field is defined by 
integrating the 3D density field along one direction, say z-axis, 

to(x) EE J p(r)dz , LOq EE p L . (22) 

The over-density in 2D space is: 

M*) = ^4^ = 7 /MO dz. (23) 

Eq. (6) gives the Fourier transformation of the over-density in 3D space. The Fourier trans- 
formation of the over-density in 2D space can be expressed as 

52 D (k 2D ) = \\ a ^ D (x)e iS ^ dx, (24) 

where A is the normalization surface area. Inserting Eq. (23) into Eq. (24), we get 

S2D(k 2D ) = W\\ S 3D (f)e i3 ^ dzdx 

J A J L 
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n J aJl - 



h J 



hD{k'2D, k z )S Dirac (k' 2 D, k 2D )5 Dira c(k z , 0) 

— # 

^3I)(^2D,0). 



(25) 



Eq. (25) gives an important hint that we can connect the statistical properties of 3D Fourier 
space with that of 2D Fourier space. Following Eq. (8), we get 



Consistent with the Cosmological Principle (the isotropic property of the Universe), Eq. (26) 
means that the 2D density field in Fourier space has the same statistical properties as the 
3D field. 

In Figure 1 we compare the 2D power spectrum with the 3D power spectrum measured 
from the simulations with initial spectral index 1, 0, —1, and —2. The number of grid 
points of the 2D FFT is 16384 2 , and the number of the 3D FFT is 1024 3 . For the 3D power 
spectrum we plot A 2 (/ci? ) (k = nkb, n — 1, 2 • • •) for the last output of the simulations until 
n = 512. As to the 2D case, we measured the power spectrum with the last three outputs of 
the simulations under the scaling transformation, and plot A 2 (kR ) ( k = 2irn/L,n = 1,2 
• • •) until n = 3000. Here we do not correct for the alias effect, and find that the alias begins 
to affect the power spectrum when the wavenumber is larger than kAL = as indicated 

by the vertical lines in the picture. This picture also shows that the 2D power spectrum is 
identical to the 3D power spectrum when the alias effect is negligible. 

In Figure 2 we compare the reduced bispectrum measured in the 2D and 3D density 
fields with the simulation of n — —1. We express the reduced bispectrum as a function of 
k, v and as defined in Eq. (9). Each panel shows Q(k,v,6) as a function of the angle 9 
between k\ and k 2 with different k and v. For the 3D bispectrum we only show the results of 
the last output for this simulation (a = 1). For the 2D case we get the reduced bispectrum 
with the last two outputs of the simulation (a = 0.792, 1) under the scaling transformation. 
The corresponding k values of the triangles are less than so the alias does not affect the 
measured results. The figure shows that the bispectrum measured in 2D agrees very well 
with that in 3D, as expected. 

From Figures 1 and 2, we conclude that it is reliable to use the 2D FFT instead of the 3D 
FFT to calculate the power spectrum and bispectrum on small scales. For the 2D particle 
distribution, an FFT with N m = 8192 requires computer memory about 0.3 Gbyte, even 



— * — * 

(S2D(k 1 )5 2 D(k 2 )) 
(5 2 D(ki)5 2 D(k 2 )5 2D (k 3 )) 



(Mfci,o)Mfc 2 ,o)), 
(Mfci,o)Mfc2,o)M*3,o)>. 



(26) 
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that with N m = 16384 requires memory about 1 Gbyte only. Therefore we can extend the 
bispectrum measurement into the strong nonlinear regime with little computer limitation. 

In fact, N m = 8192 is sufficient to avoid the alias effect. In Figure 3, we show the 
power spectra at different epochs measured with the same number of FFT grid points. The 
deviation of the power spectrum from the scaling expectation at k < Ical should be attributed 
to the numerical artifacts in the simulations. From the figure, we find that these artifacts 
begin to influence the power spectrum at k w 0.1k v for late outputs where k v = 2-k/t]. 
We believe this is mainly caused by the force softening. For early outputs, the deviation 
from the scaling happens at a larger scale, which we attribute to the cutoff of the initial 
fluctuation aX k > k Ny as well as the initial distribution of the particles on grid. All this 
shows that N m = 8192 is sufficient to explore the non-linear features in our simulations. 
In Figure 4 we check the softening effect on the reduced bispectrum. For simplicity we 
give the results only for spectral index —1. The left panels show the reduced bispectrum 
for equilateral triangles at four epochs up to the wavenumber where the softening begins to 
affect the power spectrum (as defined in Figure 3). The points with error bars are measured 
from the 2D density field, and those without error bars are obtained from the 3D density 
field. The error bars of the 2D bispectra are estimated from three projections. We compare 
these results on the right panel using the similarity scaling, which shows that the reduced 
bispectrum is not affected by the softening for the wavenumbers less than the softening limit 
set by the power spectrum analysis (Figure 3). We will use this criteria to minimize the 
softening effect in our Q measurement. 

We should point out that for a fixed range of k, the number of the Fourier modes is 
much smaller in 2D than in 3D. This limitation implies that the 2D FFT is not appropriate 
for the measurement at small wavenumbers. In our current analysis, we do both the 3D FFT 
with N m = 1024 and the 2D FFT with N m = 8192 for each simulation. For the wavenumber 

— * — * 

k < 2ir/L x 512/3, we use 5(k) based on the 3D FFT; otherwise we use S(k) based on the 
2D FFT. For the 2D case, we have projected the density field along three axes. The 2D 
bispectrum is measured by averaging over the results of the three projections. 

6. Results 

6.1. Numerical results 

In this section, we show the reduced bispectrum from the quasilinear regime to the strong 
non-linear regime measured with our new method. We have measured the bispectrum for 
many triangle shapes at different scales. In each case, we scale the results of two or three 
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outputs with the characteristic scale defined by Eq. (1) in order to make sure that these 
results are not affected by numerical effects. For simplicity, in Figure 5 we show only a few 
measured bispectra for spectral index 1, 0, —1, —2, and these results scale very well. From 
the measured results (Figures 7, 9 and 11; left columns), we find that the results are in 
good agreement with the one-loop PT prediction in the quasilinear regime (Scoccimarro et 
al. 1998): the reduced bispectrum Q is higher than the second-order PT for n = — 2 and is 
smaller for n > — 1. We also find that Q depends on the triangle shape and size even in the 
strong nonlinear regime, which is in contrast with the hypothesis adopted by SF99 that Q 
is a constant in this regime. It decreases with the initial spectral index and increases with 
the wavenumber, similar to what found for cold dark matter models (Jing & Borner 1998, 
Ma & Fry 2000, Scoccimarro, et al. 2001b, Takada & Jain 2003a). 



6.2. A fitting formula for the bispectrum 

SF99 presented a fitting formula for the bispectrum in scale-free clustering models. From 
Figures 7, 9 and 11 we find that the reduced bispectrum predicted by this formula agrees with 
the measured bispectrum at linear and quasilinear scales, and works well for some triangles 
with special shapes (such as ki = 2ki in the n = — 1 model) in the strong nonlinear regime. 
But this formula cannot generally follow the bispectrum accurately at strongly nonlinear 
scales, because they assumed that the normalized bispectrum Q is a constant for a given 
initial spectral index, independent of the triangle's wavenumber and shape. We provide a 
new fitting formula in this section to describe the nonlinear evolution of the bispectrum. 

Hamilton et al. (1991) proposed a universal empirical relation between the linear £^ and 
non- linear £tvl two-point correlation functions. The relation is a powerful tool for predicting 
£nl in cosmological models. Later, Jain, Mo, & White (1995) demonstrated that the relation 
of Hamilton et al. fails for the scale-free model of n = —2. Peacock & Dodds (1996) examined 
a large set of scale-free models and CDM models, and obtained a accurate fitting formula for 
the power spectrum which agrees with the results of N-body simulations. Motivated by their 
work, we find that the ratio of the measured bispectrum to the weakly non-linear bispectrum 
(the second-order PT) has interesting properties, especially the behavior of {Q n i/Qi) 1+0 ' 25n 
(—2 < n < 0). A few typical examples of (Q n i/Qi) 1+0 ' 25n are shown in Figures 6, 8 and 
10 (the open symbols) for the three models with n — 0, — 1, and —2. From these results 
we expect that the relation between the weakly non-linear bispectrum and the nonlinear 
bispectrum can be expressed as: 

Qni(kR , v, 9) = ti{ {1+0 - 25n) (kR , v, 9, n)Q z (^o, v,9)), (27) 

where n satisfies — 2 < n < 0. At linear scales Q n i ~ Qi, and /„; takes an asymptotic form 
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as f n i(kR 0: v,6,n) ~ 1 in the linear regime. From the measured (Q n i/Qi) 1+0 ' 25n , we can find 
that (Qni/Qi) 1+0 ' 25n as a function of 9 is approximated by a Gaussian function. Its difference 
from the Gaussian function is dependent on the triangle shape, scale and spectral index n. 
Taking into account of these factors, we propose the following form for f n i(kR ,v,9,n): 

tun ^ ( (111 -a2{kR ,v,n)) 2 \ , , 

ai(kRQ,v,n) exp — — + a 4 (kR , v, n) 

f (up „ a n \ V a 3 (kR ,v,n) J 

J n i{KK , v, V, n - — — — — — , (^8) 

1 + a 5 (kR ,v,n)(9/TT) 2 + a 6 (fci? , n)(6>/7r) 4 

where ai(kR ,v,n) are: 

ai (kR ,v,n) = {1 2 rr^—j + [(4 + n)0.002t; + 

exp(0.5(l — n)kR ) + 1 

(0.012 + 0.008n)]A; J Ro}[(0.1 - 0.3n)w + 0.4], 
a2(kRo,v,n) = 0.5 + 0.2?;, 

a 3 (kR ,v,n) = [1. 2 , + (0.01m; + 0.001)A;i? ]0.04 + 

exp(0.5/cit ) + 1 

0.06t; + 0.1 (1 - n) , 

04(^0, v,n) = {1.1 + [(n 2 ) L3 0.05 - 0.2] tanh(2A; J R ) + 0.15 exp(-(0.3A; J Ro) 2 ) + 
(Q Ql _ 0_05_OM5v )kRo}/[1 + (Q 4 _ Q 2v)1+kRo]f 

a 6 (kR ,v,n) = [0.8 + (0.01 - 0.005/ (u + 0.1))fci^] 

exp ((0.5 — n)kRo) + 1 

[2 ^— ] + 0.2ntanh(2A; J Ro - 1), 

v + 0.1 

a 5 (kR ,v,n) = —0.7a 6 (kRo,v,n), (29) 

where tanh(x) = (exp(x) — exp(— x))/(exp(x) + exp(— x)) is the hyperbolic function. 

We have obtained the best-fitting parameters in ai(kRQ, v, n) by doing a x 2 minimization 
between the predicted and measured bispectrum for all triangles. In Figures 6, 8 and 10 we 
compare the best fitting function f n i with our measured (Q n i/Qi) 1+0 ' 25n - The figures show 
that the fitting formula (28) works very well for — 2 < n < 0, for all triangle shapes, and for 
all scales (characterized by kR ) studied. With the function j nh we can convert the weakly 
non-linear bispectrum into the non-linear bispectrum. In Figures 7, 9 and 11 we compare 
our measured reduced bispectra with the prediction of our fitting formula. There we also 
plot the prediction of the fitting formula of SF99. From the figures, one can see that the 
fitting formula obtained in this paper can accurately match the simulation results, much 
better than the formula of SF99. 



To quantify the accuracy of our fitting formula, we plot the percentage of the deviation 
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AQ/Q in Figures 12, 13 and 14. The deviation is defined as: 

simu | 



(30) 



Q Q simu 

where Qfu is the prediction of our fitting formula, and Q S imu is the reduced bispectrum 
measured from the simulations. Similarly we also estimate the accuracy of the fitting formula 
of SF99 that is also shown in the figures. The figures show that our fitting formula can match 
the measured results typically at an accuracy of 10 percent. The deviation is slightly larger, 
about 20 ~ 30%, when 9 ~ it and kRo < 1, which could be attributed to some stochastic 
fluctuation in Q S i mu for the limited number of k-triangles in the configuration. We also 
see that the deviation for the SF99 formula is much larger, almost about 50 percent in the 
strongly nonlinear regime. 

Although we only have a single realization for these simulations, we measured the bis- 
pectrum for two or three outputs scaled according to the similarity solution. We regard each 
output as a realization, and estimate the errors of the bispectra from the different outputs. 
We also have taken the projections along different axes as independent realizations when we 
estimate errors for the 2D bispectra. In Figures 7, 9 and 11, the error bars are estimated 
with this method. 



7. Summary and discussion 

In this paper, we have accurately measured the bispectrum for four scale-free models. 
The measurement is based on a new method that can effectively eliminate the alias and 
numerical artifacts, and reliably extend the analysis into the strongly non-linear regime. The 
work also makes use of a set of state-of-the art N-body simulations of scale-free hierarchical 
models that have a significantly larger dynamical range than the previous studies. With these 
measured results, we demonstrated that the measured bispectrum depends on the shape and 
size of fc-triangle even in the strongly nonlinear regime. It increases with wavenumber and 
decreases with the spectral index. These results are consistent with that those found for the 
three-point correlation for CDM models (Jing & Borner 1998, Ma & Fry 2000, Scoccimarro 
& Couchman 2001 , Scoccimarro et al. 2001b, Takada & Jain 2003a), but are in contrast 
with the hypothesis that Q is a constant in the strongly non-linear regime (SF99). We also 
show that the fitting formula of SF99 does not describe our simulation results well, with a 
typical error about 40 percent. In the end, we present a new fitting formula for the reduced 
bispectrum that is valid for — 2 < n < with a typical error of 10 percent only. 

Our new method for measuring the bispectra is to use the property that the 2D power 
spectrum and bispectrum are identical to the 3D ones. This property can be easily proved 
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and has been tested with our N-body simulations. As the 2D FFT requires less computer 
memory than the 3D FFT, we can extend our analysis of the bispectrum into very small 
scales with little computer limitation. We also use the scaling properties of the scale-free 
models to correct for all known numerical artifacts and to identify those regimes where the 
bispectra can reliably measured. 

Although we have obtained an empirical formula of the bispectrum for scale-free models 
of the spectral index —2 < n < 0, some issues need to be addressed in future work. One 
obvious aspect is that the formula does not work for n = 1, because the tree-level PT cannot 
predict the weakly non-linear bispectrum for this model. Fortunately the slope of the power 
spectrum in CDM models, which are the most plausible theory for the structure formation in 
the Universe, is in the range —1 to —3 at the non-linear scales. Therefore, we may generalize 
our fitting formula to CDM models by considering the change of the power spectrum slope 
with scale as well as the deviation from the Einstein-de Sitter model. We will study the 
bispectrum in CDM models in a future paper. Another issue is that our fitting formula is 
purely empirical. Considering that our measured Q results are in good agreement with the 
one-loop PT prediction in the quasilinear regime (Scoccimarro et al. 1998) and that the halo 
model can successfully match the three-point correlation function in the strongly non-linear 
regime (Takada & Jain 2003a), we think it would be possible to combine these theoretical 
predictions to find a more theory-oriented fitting formula for Q. 

We thank Ue-Li Pen for useful discussion at the initial stage of this work. The work 
is supported in part by NKBRSF (G19990754), by NSFC (Nos.10125314, 10373012), and 
by the CAS-MPG exchange program. Numerical simulations presented in this paper were 
carried out at ADAC (the Astronomical Data Analysis Center) of the National Astronomical 
Observatory, Japan. 
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Table 1: The scale-free simulations of 512 3 particles 
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10 1 10 2 10 3 
kR„ 



Fig. 1. — Comparison of the power spectrum measured from the 2D and 3D density 
fields. The data are the last three outputs of the simulations with the spectral index 
n = 1,0,-1,-2. In each panel, the three short-dashed lines (middle) represent the 2D 
power spectrum measured from the last three outputs of the simulation, the solid line (bot- 
tom) is the 3D power spectrum measured from the last output of the simulation, and the 
long-dashed line (up) is the prediction by the fitting formula of Peacock & Dodds (1996). 
The value of NG represents the number of grid points adopted in their FFT respectively. We 
plot A 2 (kR ) (k = 2im/L, n = 1, 2 • • •.) until n = 512 for the 3D power spectrum, and until 
n = 3000 for the 2D power spectrum. The wavenumber k^L at the vertical line is k^ y /3 for 
3D power spectrum. Rq is the characteristic scale defined in Eq. (1). 
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Fig. 2. — Comparison of the reduced bispectrum measured from the 2D and the 3D density 
fields. The data are the last two outputs (a = 0.792, a = 1.0) of the simulation with spectral 
index n — — 1. The solid line corresponds to the 3D reduced bispectrum at the last output 
of the simulation, and the open symbols are for the 2D reduced bispectrum at the last two 
outputs of the simulation under the scaling transformation. The left panels give Q(k,v,9) 
for kR = 11.8, and the right ones for kR = 17.8. 
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Fig. 3. — The power spectrum measured by FFT with the same number of grid points 
(NG = 16384 2 ) for five epochs of the simulations, and scaled by the characteristic scale. 
These results are plotted by the lines A, B, C, D, E in each panel until k^i, and their epochs 
are shown in the picture. The numerical artifacts begin to affect the power spectrum at the 
vertical line indicated by k NB , k NC , k ND , k^E (in units of k v , r\ is the softening length), 
at which the deviation of the power spectrum from the true power spectrum is about 10 
percent. 
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Fig. 4. — The reduced bispectrum for equilateral triangles measured by FFT with iVG 2 D — 
8192 2 or NG30 = 1024 3 for four epochs of the n = —1 simulation. The results are plotted 
up to the wavenumber at which the softening begins to influence the power spectrum (the 
vertical lines in Figure 3). Each panel on the left shows the results for one epoch. The points 
with and without error bars are measured from the 2D and 3D density fields respectively. 
We compare these results on the right panel where for clarity we do not show the errors. 
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Fig. 5. — The reduced bispectrum Q(k,v,9) as a function of the angle 9 between k\ and hi 
for two or three outputs scaled by the characteristic scale. These results are measured by 3D 
FFT or 2D FFT as indicated at the bottom. The scales are also indicated. Different symbols 
(open circles, open triangles, open squares) stand for the results at different outputs. 
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Fig. 6. — The value of {Qni/Qi) 1 ^ 1 ^' 25 ^ as a function of 6 for the spectral index n — 0. 
Qn« is the reduced bispectrum measured from the simulation, and Qi is the weakly non- 
linear reduced bispectrum predicted by PT2. The solid line shows our fitting formula for 
(Qni/Qi) 1 ^ 1+0 ' 25n \ and the circle symbols are the simulation results. Vertical panel columns 
have the same kR , as indicated at the bottom. 
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Fig. 7. — The reduced bispectrum Q(k, v, 6) measured from the simulation with the spectral 
index n = (open circles), compared with the predictions by PT2 (long-dashed lines), by 
the fitting formula of SF99 (short-dashed lines), and by the fitting formula in this paper 
(solid lines). Vertical panel columns have the same kR , as indicated at the bottom. 
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Fig. 8. — Same as Fig. 6, but for the spectral index n = 



-1. 
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Fig. 9. — Same as Fig. 7, but for the spectral index n = 



-1. 
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Fig. 10. — Same as Fig.6, but for the spectral index n = 



-2. 
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Fig. 11.- 



Same as Fig. 7, but for the spectral index n = 



-2. 
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Fig. 12. — The relative accuracy of the fitting formulae for the bispectrum for the spectral 
index n — 0. The open circles are for the formula of SF99, and the solid lines are for the 
formula obtained in this paper. 
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Fig. 13.- 



Same as Fig. 12, but for the spectral index n = 



-1. 
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Fig. 14. — Same as Fig. 12, but for the spectral index n = 
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